function [R_maneuver,V_maneuver,t_maneuver,delV] = Run_StationKeeping_LT(S_chase_eci,S_target_eci,S_target_ecef,...
    t_crossing,UTC_mission_start,targetAccel,mass,isp,T);
%% Constants
options = bvpset('RelTol', 1e-6, 'AbsTol', 1e-9,'NMax',5000);
mu = 398600.4418;
mu_unit = 1;
siderial = 86164.1;
sma_geo = (mu*(siderial/(2*pi))^2)^(1/3);
We = (2*pi / siderial) * [0;0;1];
options = odeset('RelTol' , 1E-6);
%% Find the states of the chase vehicle in HCW frame of the S_target_eci;
% Define new variables
year = UTC_mission_start(1);
mon  = UTC_mission_start(2);
day  = UTC_mission_start(3);
hour = UTC_mission_start(4);
minute  = UTC_mission_start(5);
sec  = UTC_mission_start(6);

%% Run low thrust maneuver using homotopy method.

% 1. Solve tpbvp until the target accel is met

% Make initial guess vector
DU2km = sma_geo;
TU2s = sqrt(DU2km^3/mu);
VU2kms = DU2km/TU2s;
target_Accel = targetAccel*(TU2s^2/(DU2km*1000));
isp_TU = isp*(1/TU2s);
g_DUTU2 = 9.81*(1/DU2km)*(1/1000)*(TU2s^2);
cexh = isp_TU*g_DUTU2;

COE_chase           = State2Coe(S_chase_eci,mu);
S_chaseCheck   = Coe2State([COE_chase(1:5);COE_chase(6)],mu);
%%%%%

el0.a               = COE_chase(1);
el0.e               = COE_chase(2);
el0.i               = COE_chase(3)*pi/180;
el0.aop             = COE_chase(4)*pi/180;
el0.raan            = COE_chase(5)*pi/180;
checker = elements2rtp(el0,COE_chase(6)*pi/180,mu);


%%%%%


COE_chase(1)        = COE_chase(1)/DU2km;
el0.a               = COE_chase(1);
el0.e               = COE_chase(2);
el0.i               = COE_chase(3)*pi/180;
el0.aop             = COE_chase(4)*pi/180;
el0.raan            = COE_chase(5)*pi/180;

chase_ta0 = COE_chase(6)*pi/180;

COE_target      = State2Coe(S_target_eci,mu);
target_ta0      = atan2(S_target_eci(2),S_target_eci(1));
if target_ta0 < 0
    target_ta0  = target_ta0+2*pi;
end

% S_targetCheck   = Coe2State([COE_target(1:3);0;0;COE_target(6)],mu);
COE_target(1)   = COE_target(1)/DU2km;
el1.a           = COE_target(1);
el1.e           = COE_target(2);
el1.i           = 0;
el1.aop         = 0;
el1.raan        = 0;

S_o             = [elements2rtp_xyzvxyz(el0,chase_ta0,mu_unit,S_chase_eci(1:3)/DU2km,S_chase_eci(4:6)/VU2kms);target_Accel];

ft_guess = 2*pi;

[XT, XS]    = ode45(@(t, s) xeoms(t, s, mu_unit, cexh), linspace(0, ft_guess, 200), S_o);
x           = XS(:, 1).*sin(XS(:, 2)).*cos(XS(:, 3));
y           = XS(:, 1).*sin(XS(:, 2)).*sin(XS(:, 3));
z           = XS(:, 1).*cos(XS(:, 2));
target_taf  = ft_guess+target_ta0;
if target_taf > 2*pi
    
    target_taf = target_taf-2*pi;
    
end

eci2ecef = Eci2Ecef([UTC_mission_start(1:5),UTC_mission_start(6)+t_crossing+(ft_guess*TU2s)]);
Seci_target_tf = eci2ecef'*S_target_ecef(1:3);

anglecheck = atan2(Seci_target_tf(2),Seci_target_tf(1));
if anglecheck < 0
    
    anglecheck = anglecheck+2*pi;
    
end

L0 = [0 0 0 1 0 -1 1];
[~, LS] = ode45(@(t, s) leoms(t, s, mu_unit, cexh, XT, XS), XT(end:-1:1), L0);
LS = LS(end:-1:1, :);

[~, XS0] = ode45(@(t, s) xeoms(t, s, mu_unit, cexh), linspace(0, ft_guess, 101), [elements2rtp(el0, chase_ta0, mu_unit)' 0]);
x0 = XS0(:, 1).*sin(XS0(:, 2)).*cos(XS0(:, 3));
y0 = XS0(:, 1).*sin(XS0(:, 2)).*sin(XS0(:, 3));
z0 = XS0(:, 1).*cos(XS0(:, 2));

[~, XS1] = ode45(@(t, s) xeoms(t, s, mu_unit, cexh), linspace(0, ft_guess, 101), [elements2rtp(el1, target_ta0, mu_unit)' 0]);
x1 = XS1(:, 1).*sin(XS1(:, 2)).*cos(XS1(:, 3));
y1 = XS1(:, 1).*sin(XS1(:, 2)).*sin(XS1(:, 3));
z1 = XS1(:, 1).*cos(XS1(:, 2));
% Plot things
figure;
axis equal
hold on
plot3(x, y, z,'b',...
    x(1), y(1), z(1),'bo',...
    x(end), y(end), z(end),'bo');
plot3(x0, y0, z0,'r',...
    x0(1), y0(1), z0(1),'ro',...
    x0(end), y0(end), z0(end),'ro');
plot3(x1, y1, z1,'k',...
    x1(1), y1(1), z1(1),'ko',...
    x1(end), y1(end), z1(end),'ko');
solinit.x = XT;
solinit.y = [XS LS]';
sol = bvp4c(@(t, s) eleoms(t, s, mu_unit, cexh), ...
    @(sa, sb) bcfunc_xyz(sa, sb, el0, el1, mu_unit,chase_ta0,S_chase_eci(1:3)/DU2km,target_Accel,S_chase_eci(4:6)/VU2kms),...
    solinit, ...
    options);
target_ta = sol.y(7,1);
r = sol.y(1, :);
theta = sol.y(2, :);
phi = sol.y(3, :);
x = r.*sin(theta).*cos(phi);
y = r.*sin(theta).*sin(phi);
z = r.*cos(theta);
plot3(x, y, z,'c','LineWidth',3)
plot3(x(1), y(1), z(1),'co',...
    x(end), y(end), z(end),'cd');
if (sol.y(3,end)-target_taf) > pi
    target_taf = target_taf+ft_guess;
end


solinit.x = sol.x;
solinit.y = sol.y;
sol = bvp5c(@(t, s) eleoms(t, s, mu_unit, cexh),...
    @(sa, sb) bcfunc_SK_xyz(sa, sb, el0, el1, mu_unit,target_taf,chase_ta0,S_chase_eci(1:3)/DU2km,target_ta,S_chase_eci(4:6)/VU2kms),...
    solinit, options);
while sol.stats.maxerr > 5e-5
    solinit.x = sol.x;
    solinit.y = sol.y;
    sol = bvp5c(@(t, s) eleoms(t, s, mu_unit, cexh),...
        @(sa, sb) bcfunc_SK_xyz(sa, sb, el0, el1, mu_unit,target_taf,chase_ta0,S_chase_eci(1:3)/DU2km,target_ta,S_chase_eci(4:6)/VU2kms),...
        solinit, options);
end
% anglediff = phi_f - target_taf;
r = sol.y(1, :);
theta = sol.y(2, :);
phi = sol.y(3, :);
x = r.*sin(theta).*cos(phi);
y = r.*sin(theta).*sin(phi);
z = r.*cos(theta);
plot3(x, y, z,'m','LineWidth',3)
plot3(x(1), y(1), z(1),'mo',...
    x(end), y(end), z(end),'md');   

% 1. Change the state in ECI frame.

r = sol.y(1, :);
theta = sol.y(2, :);
phi = sol.y(3, :);

x = r.*sin(theta).*cos(phi);
y = r.*sin(theta).*sin(phi);
z = r.*cos(theta);

R = [x;y;z]*DU2km;

for ii = 1:length(R)
    theta = sol.y(2, ii);
    phi = sol.y(3, ii);
    V_rtp = sol.y(4:6,ii);
    Axyz2rtp = [sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta);
        cos(theta)*cos(phi), cos(theta)*sin(phi), -sin(theta);
        -sin(phi), cos(phi), 0];
    V(:,ii) = (Axyz2rtp'*V_rtp)*VU2kms;
end



t = t_crossing+0.000001 + sol.x*TU2s;

Rf_sol = R(:,end);
eci2ecef = Eci2Ecef([UTC_mission_start(1:5),UTC_mission_start(6)+t(end)]);
Rf_target = eci2ecef'*S_target_ecef(1:3);
posdiff = norm(Rf_sol-Rf_target);



%% Save history
R_maneuver = [R];
% V_maneuver = [V_eci(:,1:n-1),Vt(:,end)+0.01*rand];
V_maneuver = [V];

t_maneuver = [t];

ta_dimension = sol.y(7,:)*(1/TU2s)^2*DU2km*1000;
delV = abs(trapz(t,ta_dimension));




end

